###################################################################################################
##                                                                                               ##
##  This is the replication file for the analysis in Tsai and Lin. 2017. "Modeling Guessing      ##
##  Components in the Measurement of Political Knowledge."                                       ##
##                                                                                               ##
###################################################################################################
########################    APPENDIX A    ##############################################
##
##    Please see the R file "ANES_eggs4_replication"
##
########################    APPENDIX B    ##############################################
########################    CONVERGENCE TESTS FOR MCMC    ##############################################
## LOAD PACKAGES
library(R2jags);
library(scales);

## LOAD REQUIRED FUNCTIONS
source("/knowledge_functions.R");

## LOAD THE DATA
mod1 <- read.csv(file="2pl_para_mul.csv", header=TRUE);
dim(mod1);

mod2 <- read.csv(file="3pl_para_mul.csv", header=TRUE);
dim(mod2);

mod3 <- read.csv(file="guess_para_mul.csv", header=TRUE);
dim(mod3);


## transform to MCMC object
chain.1 <- as.mcmc(mod3[1:5000,]);
chain.2 <- as.mcmc(mod3[5001:10000,]);
chain.3 <- as.mcmc(mod3[10001:15000,]);
dim(chain.1);

chains <- as.mcmc.list(list(chain.1, chain.2, chain.3));

########################  TABLE B-1: GELMAN-RUBIN TEST  #######################
gelman.diag(chains);

########################  TABLE B-2: GEWEKE TEST  #######################
geweke.diag(chains);


###################  PARAMETERS OF PLOTS
par(mfrow=c(1,1),mar=c(4,5,1,1),oma=c(1,1,1,1));

########################  FIGURE B-1a: TRACE PLOT OF b[1]  #######################
plot(mod3[1:5000,"b.1."], type="n", xlab="Number of Iterations", ylab=expression(paste("Estimate of ",b[1])), ylim=c(0,1), yaxt="n", xaxt="n", bty="n", cex.lab=1.2);

axis(side=1, las=1, tick=T, cex.axis=1.1, at=seq(0,5000,1000) ); # X-axis
axis(side=2, las=1, tick=T, cex.axis=1.1, at=seq(0,1,0.2) ); # Y-axis

lines(mod3[5001:10000,"b.1."], col=alpha("blue", 0.4) );
lines(mod3[10001:15000,"b.1."], col=alpha("red", 0.4) );
lines(mod3[1:5000,"b.1."], col=alpha("black", 0.8) );


########################  FIGURE B-1b: TRACE PLOT OF b[2]  #######################
plot(mod3[1:5000,"b.2."], type="n", xlab="Number of Iterations", ylab=expression(paste("Estimate of ",b[2])), ylim=c(0,1), yaxt="n", xaxt="n", bty="n", cex.lab=1.2);

axis(side=1, las=1, tick=T, cex.axis=1.1, at=seq(0,5000,1000) ); # X-axis
axis(side=2, las=1, tick=T, cex.axis=1.1, at=seq(0,1,0.2) ); # Y-axis

lines(mod3[5001:10000,"b.2."], col=alpha("blue", 0.4) );
lines(mod3[10001:15000,"b.2."], col=alpha("red", 0.4) );
lines(mod3[1:5000,"b.2."], col=alpha("black", 0.8) );


########################  FIGURE B-1c: TRACE PLOT OF b[3]  #######################
plot(mod3[1:5000,"b.3."], type="n", xlab="Number of Iterations", ylab=expression(paste("Estimate of ",b[3])), ylim=c(0,1), yaxt="n", xaxt="n", bty="n", cex.lab=1.2);

axis(side=1, las=1, tick=T, cex.axis=1.1, at=seq(0,5000,1000) ); # X-axis
axis(side=2, las=1, tick=T, cex.axis=1.1, at=seq(0,1,0.2) ); # Y-axis

lines(mod3[5001:10000,"b.3."], col=alpha("blue", 0.4) );
lines(mod3[10001:15000,"b.3."], col=alpha("red", 0.4) );
lines(mod3[1:5000,"b.3."], col=alpha("black", 0.8) );


########################  FIGURE B-1d: TRACE PLOT OF b[4]  #######################
plot(mod3[1:5000,"b.4."], type="n", xlab="Number of Iterations", ylab=expression(paste("Estimate of ",b[4])), ylim=c(0,1), yaxt="n", xaxt="n", bty="n", cex.lab=1.2);

axis(side=1, las=1, tick=T, cex.axis=1.1, at=seq(0,5000,1000) ); # X-axis
axis(side=2, las=1, tick=T, cex.axis=1.1, at=seq(0,1,0.2) ); # Y-axis

lines(mod3[5001:10000,"b.4."], col=alpha("blue", 0.4) );
lines(mod3[10001:15000,"b.4."], col=alpha("red", 0.4) );
lines(mod3[1:5000,"b.4."], col=alpha("black", 0.8) );


##
mcmc.mod3 <- as.mcmc(mod3);
traceplot(mcmc.mod3[,10]);

###################################################################################################
########################    APPENDIX C    ##############################################
######################################################
##
##    Please see line 583 in the R file "main_analysis"
##
###################################################################################################
########################    APPENDIX D    ##############################################
######################################################
## LOAD THE DATA
theta.2pl <- read.csv(file="2pl_theta_mul.csv", header=TRUE);
dim(theta.2pl);
mu.2pl <- apply(theta.2pl, 2, mean);

theta.3pl <- read.csv(file="3pl_theta_mul.csv", header=TRUE);
dim(theta.3pl);
mu.3pl <- apply(theta.3pl, 2, mean);

theta.g <- read.csv(file="guess_theta_mul.csv", header=TRUE);
dim(theta.g);
mu.g <- apply(theta.g, 2, mean);

sort.index <- names(sort(mu.g));


###################  PARAMETERS OF PLOTS
par(mfrow=c(1,1),mar=c(4,4,1,1),oma=c(1,1,1,1));


###################  FIGURE D1-1: theta.guess against theta.2pl
plot(mu.g[sort.index], mu.2pl[sort.index], type="n", xlab="Knowledge Estimates from 2PL-G", ylab="Knowledge Estimates from 2PL", main="Knowledge Estimates from 2PL-G and 2PL", xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), yaxt="n", xaxt="n", bty="n");

axis(side=2, las=1, tick=FALSE, cex.axis=1.2, at=seq(-2.5,2.5,0.5) ); # Y-axis
axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(-2.5,2.5,0.5) ); # X-axis

abline(v=seq(-6,6,0.01), col="gray97", lty=1, lwd=1); # background

abline(h=seq(-5,5,1), col="white", lty=1, lwd=2);
abline(v=seq(-5,5,1), col="white", lty=1, lwd=2);

abline(a=0, b=1, col="gray70", lty=2);

points(mu.g[sort.index], mu.2pl[sort.index], pch=19, col=adjustcolor("gray20", alpha.f=0.6), cex=1.25);


###################  FIGURE D1-2: theta.guess against theta.3pl
plot(mu.g[sort.index], mu.3pl[sort.index], type="n", xlab="Knowledge Estimates from 2PL-G", ylab="Knowledge Estimates from 3PL", main="Knowledge Estimates from 2PL-G and 3PL", xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5), yaxt="n", xaxt="n", bty="n");

axis(side=2, las=1, tick=FALSE, cex.axis=1.2, at=seq(-2.5,2.5,0.5) ); # Y-axis
axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(-2.5,2.5,0.5) ); # X-axis

abline(v=seq(-6,6,0.01), col="gray97", lty=1, lwd=1); # background

abline(h=seq(-5,5,1), col="white", lty=1, lwd=2);
abline(v=seq(-5,5,1), col="white", lty=1, lwd=2);

abline(a=0, b=1, col="gray70", lty=2);

points(mu.g[sort.index], mu.3pl[sort.index], pch=19, col=adjustcolor("gray20", alpha.f=0.6), cex=1.25);


###################################################################################################
########################    APPENDIX E    ##############################################
######################################################
##
##  LOAD AND REOPORT THE RESULTS
##
######################################################
## LODING FUNCTIONS
source("/knowledge_functions.R");

## LOAD THE DATA
mod1 <- read.csv(file="2pl_para2.csv", header=TRUE);
dim(mod1);

mod2 <- read.csv(file="3pl_para2.csv", header=TRUE);
dim(mod2);

mod3 <- read.csv(file="guess_para2.csv", header=TRUE);
dim(mod3);



###################  FIGURE E1: densities of ideal point estimates
###################  PARAMETERS OF PLOTS
par(mfrow=c(1,1), mar=c(4,4,3,1), oma=c(1,1,1,1));

plot(mu.g, dnorm(mu.g), type="n", xlab="Political Knowledge", ylab="Density", main="Distributions of Political Knowledge", xlim=c(-2,2), ylim=c(0,1), yaxt="n", xaxt="n", bty="n");

axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(-2,2,0.5) ); # X-axis
axis(side=2, las=1, tick=FALSE, cex.axis=1.2, at=seq(0,1,0.5) ); # Y-axis

abline(v=seq(-2.5,2.5,0.01), col="gray97", lty=1, lwd=1);

abline(h=c(0,0.5,1), col="white", lty=1, lwd=3);

abline(v=seq(-2,2,0.5), col="gray100", lty=1, lwd=1);
abline(v=seq(-2,2,1), col="white", lty=1, lwd=3);

# 2PL-G
lines(density(mu.g, adjust=2), col="black", lty=1, lwd=2, type="l");
#lines(density(mu.g, adjust=2), col=alpha("black", 0.2), lty=1, lwd=2, type="h");

# 2PL
lines(density(mu.2pl, adjust=2), col="black", lty=2, lwd=2, type="l");
#lines(density(mu.2pl, adjust=2), col=alpha("red", 0.1), lty=1, lwd=2, type="h");

# 3PL
lines(density(mu.3pl, adjust=2), col="black", lty=3, lwd=2, type="l");
#lines(density(mu.3pl, adjust=2), col=alpha("blue", 0.1), lty=1, lwd=2, type="h");

#
legend(x=0.5, y=0.75, legend=c("2PL-G", "3PL", "2PL"), lty=c(1,3,2), cex=1, lwd=2, bty="n" );


#######################################################################################
######################    Comparisons of Models    ######################
## LOAD THE DATA
vp_2012 <- read.csv(file="/teds2012_ind.csv", header=TRUE);
dim(vp_2012);

res.data <- vp_2012[,c("re_G01","re_G02","re_G03","re_G04","re_G05","re_G06","re_G07","G04","G05","G06","G07")];

attach(vp_2012);

res.matrix <- cbind(re_G06, re_G04, re_G05, re_G07);
dim(res.matrix);

detach(vp_2012);

##  Distribution of Choice Options
res.data$n.correct <- rowSums(res.data[,c(4:7)], na.rm=TRUE);

index <- which(res.data$n.correct >= 1 & res.data$n.correct <= 3);


##  2PL: The Proportion of Correct Predictions
n.correct <- data.frame(matrix(NA, nrow=dim(mod1)[1], ncol=8));
colnames(n.correct) <- c("m.q6", "m.q4", "m.q5", "m.q7", "n.q6", "n.q4", "n.q5", "n.q7");

for (i in 1:dim(mod1)[1]) {
	pred.2pl <- pred.rate(response=res.matrix[index,], theta=theta.2pl[i, index], beta=mod1[i,6:9], alpha=mod1[i,1:4], model="2PL");
	
	n.correct[i,1:4] <- pred.2pl$item.match;
	n.correct[i,5:8] <- pred.2pl$n.answer;
}

##  Correct Predictions by Items
n.correct[,9:12] <- n.correct[,1:4]/n.correct[,5:8];

par.mod1 <- cbind(round(apply(n.correct[,9:12], 2, mean), 3), HPDint(n.correct[,9:12], prob=0.90)$HPDinterval);
par.mod1;


##  3PL: The Proportion of Correct Predictions
n.correct2 <- data.frame(matrix(NA, nrow=dim(mod1)[1], ncol=8));
colnames(n.correct2) <- c("m.q6", "m.q4", "m.q5", "m.q7", "n.q6", "n.q4", "n.q5", "n.q7");

for (i in 1:dim(mod1)[1]) {
	pred.3pl <- pred.rate(response=res.matrix[index,], theta=theta.3pl[i, index], beta=mod2[i,6:9], alpha=mod2[i,1:4], model="3PL", c=mod2[i,10:13]);
	
	n.correct2[i,1:4] <- pred.3pl$item.match;
	n.correct2[i,5:8] <- pred.3pl$n.answer;
}

##  Correct Predictions by Items
n.correct2[,9:12] <- n.correct2[,1:4]/n.correct2[,5:8];

par.mod2 <- cbind(round(apply(n.correct2[,9:12], 2, mean), 3), HPDint(n.correct2[,9:12], prob=0.90)$HPDinterval);
par.mod2;


##  2PL-G: The Proportion of Correct Predictions
n.correct3 <- data.frame(matrix(NA, nrow=dim(mod1)[1], ncol=8));
colnames(n.correct3) <- c("m.q6", "m.q4", "m.q5", "m.q7", "n.q6", "n.q4", "n.q5", "n.q7");

for (i in 1:dim(mod1)[1]) {
	pred.g <- pred.rate(response=res.matrix[index,], theta=theta.g[i, index], beta=mod3[i,6:9], alpha=mod3[i,1:4], model="2PL-G", b=mod3[i,10:13], M=4);
	
	n.correct3[i,1:4] <- pred.g$item.match;
	n.correct3[i,5:8] <- pred.g$n.answer;
}

##  Correct Predictions by Items
n.correct3[,9:12] <- n.correct3[,1:4]/n.correct3[,5:8];

par.mod3 <- cbind(round(apply(n.correct3[,9:12], 2, mean), 3), HPDint(n.correct3[,9:12], prob=0.90)$HPDinterval);
par.mod3;


########################  PLOT ESTIMATES  #######################
var.label <- c("Second Party", "Fin. Minister", "Unemployment", "UN Secretary");

##
par(mfrow=c(1,1), mar=c(2,4,5,2), oma=c(2,2,0,0));

########################  FIGURE E2:   #######################
plot(1:17 ~ 0, type="n", xlab="", ylab="", main="Proportions of Correct Prediction", axes=FALSE, xlim=c(0.49,1.01), cex.main=1.1);

axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(0,1,0.1) ); # X-axis
axis(side=3, las=1, tick=FALSE, cex.axis=1.2, at=seq(0,1,0.1) ); # X-axis

axis(side=2, las=1, tick=FALSE, cex.axis=0.8, at=c(3,7,11,15), label=rev(var.label) ); # Y-axis

abline(v=seq(-1,1.05,0.001), col="gray97", lty=1, lwd=1); # background

abline(h=c(2:4,6:8,10:12,14:16), lty=2, lwd=1, col="white");

abline(v=seq(0,1,0.1), col="white", lty=1, lwd=2);
abline(v=seq(0,1,0.05), col="white", lty=1, lwd=1);


# 90% credible intervals
segments(par.mod1[1:4,2], c(16,12,8,4), par.mod1[1:4,3], c(16,12,8,4), lty=6, col="gray20", lwd=3);
segments(par.mod2[1:4,2], c(15,11,7,3), par.mod2[1:4,3], c(15,11,7,3), lty=1, col="black", lwd=3);
segments(par.mod3[1:4,2], c(14,10,6,2), par.mod3[1:4,3], c(14,10,6,2), lty=2, col="gray70", lwd=3);

# Means
points(x=par.mod1[1:4,1], y=c(16,12,8,4), pch=23, cex=1.2, col="gray20", bg="gray20");
points(x=par.mod2[1:4,1], y=c(15,11,7,3), pch=22, cex=1.2, col="black", bg="black");
points(x=par.mod3[1:4,1], y=c(14,10,6,2), pch=21, cex=1.2, col="gray70", bg="gray70");

legend(x=0.5, y=17, legend=c("2PL", "3PL", "2PL-G"), lty=c(6,1,2), pch=c(18,15,16), col=c("gray20","black","gray70"), cex=1, horiz=FALSE, lwd=2, bg="gray97", bty="n" );


###################################################################################################
########################    APPENDIX F    ##############################################
#################################################################################
## LOAD THE DATA
vp_2012 <- read.csv(file="/teds2012_ind.csv", header=TRUE);
dim(vp_2012);

##
subject.info <- data.frame(id=vp_2012$id, female=vp_2012$female, party.id3=vp_2012$party.id3, edu5=vp_2012$re_edu);

female.index <- which(subject.info$female == 1);
male.index <- which(subject.info$female == 0);


##  FEMALE: Sum of the Number of Correct Responses Excluding The Given Item
res.data <- vp_2012[female.index,c("re_G01","re_G02","re_G03","re_G04","re_G05","re_G06","re_G07","G04","G05","G06","G07")];

##  Distribution of Choice Options: Fin. Minister
res.data$n.correct.rm4 <- rowSums(res.data[,c(1:3,5:7)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm4", col.vars="G04", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("Jiang Yi-huah", "Chen Chun", "Mao Chi-kuo", "Lee Sush-der", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


# Unemployment
res.data$n.correct.rm5 <- rowSums(res.data[,c(1:4,6:7)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm5", col.vars="G05", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("2.3%", "4.3%", "6.3%", "8.3%", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


# Second Party
res.data$n.correct.rm6 <- rowSums(res.data[,c(1:5,7)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm6", col.vars="G06", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("KMT", "DPP", "PFP", "Non-Partisan Solidarity Union", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


# UN secretary
res.data$n.correct.rm7 <- rowSums(res.data[,c(1:6)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm7", col.vars="G07", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("Kofi Annan", "Kurt Waldheim", "Ban Ki-Moon", "Boutros Boutros-Ghali", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


##  MALE: Sum of the Number of Correct Responses Excluding The Given Item
res.data <- vp_2012[male.index,c("re_G01","re_G02","re_G03","re_G04","re_G05","re_G06","re_G07","G04","G05","G06","G07")];


##  Distribution of Choice Options: Fin. Minister
res.data$n.correct.rm4 <- rowSums(res.data[,c(1:3,5:7)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm4", col.vars="G04", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("Jiang Yi-huah", "Chen Chun", "Mao Chi-kuo", "Lee Sush-der", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


# Unemployment
res.data$n.correct.rm5 <- rowSums(res.data[,c(1:4,6:7)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm5", col.vars="G05", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("2.3%", "4.3%", "6.3%", "8.3%", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


# Second Party
res.data$n.correct.rm6 <- rowSums(res.data[,c(1:5,7)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm6", col.vars="G06", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("KMT", "DPP", "PFP", "Non-Partisan Solidarity Union", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


# UN secretary
res.data$n.correct.rm7 <- rowSums(res.data[,c(1:6)], na.rm=TRUE);
ctable <- crosstab(res.data, row.vars="n.correct.rm7", col.vars="G07", type="f");
rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );
ct <- rbind("0-1"=colSums(ctable$table[1:2,]), "2-4"=colSums(ctable$table[3:5,]), "5-6"=colSums(ctable$table[6:7,]) );

ct <- cbind(ct[,1:4], ct[,5]+ct[,6], ct[,7]);
colnames(ct) <- c("Kofi Annan", "Kurt Waldheim", "Ban Ki-Moon", "Boutros Boutros-Ghali", "Don't Know", "Sum");

ct;
round(ct/ct[,6]*100, 2);


###################################################################################################
########################    APPENDIX G    ##############################################
#################################################################################
########################    SAN MARTIN'S MODEL    #######################
## LOAD THE DATA
egg4 <- read.csv(file="/anes_specialstudies_2012egss4.csv", header=TRUE);
dim(egg4);

##
attach(egg4);

res.matrix <- cbind(re_zh1, re_zh2, re_zh3, re_zh4);
dim(res.matrix);

detach(egg4);

## SPECIFY DATA PARAMETERS
y <- res.matrix;
N <- dim(y)[1];
K <- dim(y)[2];


## SET THE SEED, WHICH TELLS US WHERE THE RANDOM PROCESS STARTS.
set.seed(07102014);

data1 <- list("y", "N", "K");

# INITIAL VALUES; LARGE DISPREAD
inits1 <- list(
	list(theta=rep(0,N), alpha=rep(0,K), tau.alpha=0.1, b=1, gamma=rep(3,K) ),
	list(theta=rep(-3,N), alpha=rep(3,K), tau.alpha=1, b=0.1, gamma=rep(0,K) ),
	list(theta=rep(3,N), alpha=rep(-3,K), tau.alpha=5, b=5, gamma=rep(-3,K) )
    );

parameters1 <- c("theta", "alpha", "sigma.alpha", "b", "gamma");

##  MARKOV CHAIN MONTE CARLO
n.chains <- length(inits1);
n.iteration <- 100000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 10;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin;  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains;  # TOTAL SAMPLES SAVED

## LOAD PACKAGES
library(R2jags);

## HAND-OFF TO JAGS
bml1 <- jags(model.file="sanmartin.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

print(bml1);


##  PROCESS THE MCMC RESULTS
mcmc.out <- as.mcmc(bml1);


##  SAVE THE RESULTS
out.matrix <- as.matrix(mcmc.out);
dim(out.matrix);
#out.matrix[1,1:50];


alpha <- out.matrix[,c(1:4,11)];
b.gamma <- out.matrix[,c(5,7:10)];
theta <- out.matrix[,12:1325];

out.para <- cbind(alpha, b.gamma);
dim(out.para);


##  Order "theta" by 1,2,3,... rather than 1,10,100,...
step1 <- sapply(strsplit(x=colnames(theta), split='[', fixed=TRUE), function(x) (x[2]));
step2 <- sapply(strsplit(x=step1, split=']', fixed=TRUE), function(x) (x[1]));
step3 <- data.frame(index=1:length(colnames(theta)), theta=as.numeric(step2));
index <- step3[order(step3[,2]), ][,1];

theta.sort <- theta[,index];


##  SAVE THE RESULTS
#write.csv(out.para, "/sanmartin_para_anes.csv", row.names=FALSE);
#write.csv(theta.sort, "/sanmartin_theta_anes.csv", row.names=FALSE);
#########################################################################

######################################################
##
##  LOAD AND REOPORT THE RESULTS
##
######################################################
## LODING FUNCTIONS
source("/knowledge_functions.R");

## LOAD THE DATA
mod1 <- read.csv(file="sanmartin_para.csv", header=TRUE);
dim(mod1);


########################  PLOT ESTIMATES  #######################
par.mod2 <- cbind(round(apply(mod1, 2, mean), 3), HPDint(mod1, prob=0.90)$HPDinterval);

var.label <- c("Fin. Minister", "Unemployment", "Second Party", "UN Secretary");

##
par(mfrow=c(1,1), mar=c(2,4,5,2), oma=c(2,2,0,0));

########################  FIGURE G1-1: ITEM DIFFICULTY FOR ITEMS 4-7  #######################
plot(1:17 ~ 0, type="n", xlab="", ylab="", main="Item Difficulty", axes=FALSE, xlim=c(-5,2), cex.main=1.1);

axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(-5,2,1) ); # X-axis
axis(side=3, las=1, tick=FALSE, cex.axis=1.2, at=seq(-5,2,1) ); # X-axis

axis(side=2, las=1, tick=FALSE, cex.axis=0.8, at=c(3,7,11,15), label=rev(var.label) ); # Y-axis

abline(v=seq(-6,6,0.01), col="gray97", lty=1, lwd=1); # background

abline(h=c(2:4,6:8,10:12,14:16), lty=2, lwd=0.5, col="white");

abline(v=seq(-5,5,1), col="white", lty=1, lwd=2);


# 90% credible intervals
segments(par.mod2[1:4,2], c(15,11,7,3), par.mod2[1:4,3], c(15,11,7,3), lty=1, col="black", lwd=3);

# Means
points(x=par.mod2[1:4,1], y=c(15,11,7,3), pch=22, cex=1.2, col="black", bg="black");


########################  FIGURE G1-2: ITEM DISCRIMINATION FOR ITEMS 4-7  #######################
plot(1:17 ~ 0, type="n", xlab="", ylab="", main="Guessing Parameter", axes=FALSE, xlim=c(-20,4), cex.main=1.1);

axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(-20,5,5) ); # X-axis
axis(side=3, las=1, tick=FALSE, cex.axis=1.2, at=seq(-20,5,5) ); # X-axis

axis(side=2, las=1, tick=FALSE, cex.axis=0.8, at=c(3,7,11,15), label=rev(var.label) ); # Y-axis

abline(v=seq(-20,6,0.01), col="gray97", lty=1, lwd=1);

abline(h=c(2:4,6:8,10:12,14:16), lty=2, lwd=0.5, col="white");

#abline(v=seq(-5,5,0.5), col="white", lty=1, lwd=1);
abline(v=seq(-20,5,5), col="white", lty=1, lwd=2);


# 90% credible intervals
segments(par.mod2[7:10,2], c(15,11,7,3), par.mod2[7:10,3], c(15,11,7,3), lty=1, col="black", lwd=3);

# Means
points(x=par.mod2[7:10,1], y=c(15,11,7,3), pch=22, cex=1.2, col="black", bg="black");



################################################################################################
## LOAD THE DATA
mod1 <- read.csv(file="sanmartin_para_anes.csv", header=TRUE);
dim(mod1);


########################  PLOT ESTIMATES  #######################
par.mod2 <- cbind(round(apply(mod1, 2, mean), 3), HPDint(mod1, prob=0.90)$HPDinterval);

var.label <- c("Chief Justice", "UK Primier", "Speaker of HOR", "Least Spending");

##
par(mfrow=c(1,1), mar=c(2,4,5,2), oma=c(2,2,0,0));


########################  FIGURE G2-1: ITEM DIFFICULTY FOR ITEMS 4-7  #######################
plot(1:17 ~ 0, type="n", xlab="", ylab="", main="Item Difficulty", axes=FALSE, xlim=c(-5,2), cex.main=1.1);

axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(-5,2,1) ); # X-axis
axis(side=3, las=1, tick=FALSE, cex.axis=1.2, at=seq(-5,2,1) ); # X-axis

axis(side=2, las=1, tick=FALSE, cex.axis=0.8, at=c(3,7,11,15), label=rev(var.label) ); # Y-axis

abline(v=seq(-6,6,0.01), col="gray97", lty=1, lwd=1); # background

abline(h=c(2:4,6:8,10:12,14:16), lty=2, lwd=0.5, col="white");

abline(v=seq(-5,5,1), col="white", lty=1, lwd=2);


# 90% credible intervals
segments(par.mod2[1:4,2], c(15,11,7,3), par.mod2[1:4,3], c(15,11,7,3), lty=1, col="black", lwd=3);

# Means
points(x=par.mod2[1:4,1], y=c(15,11,7,3), pch=22, cex=1.2, col="black", bg="black");


########################  FIGURE G2-2: ITEM DISCRIMINATION FOR ITEMS 4-7  #######################
plot(1:17 ~ 0, type="n", xlab="", ylab="", main="Guessing Parameter", axes=FALSE, xlim=c(-25,4), cex.main=1.1);

axis(side=1, las=1, tick=FALSE, cex.axis=1.2, at=seq(-25,5,5) ); # X-axis
axis(side=3, las=1, tick=FALSE, cex.axis=1.2, at=seq(-25,5,5) ); # X-axis

axis(side=2, las=1, tick=FALSE, cex.axis=0.8, at=c(3,7,11,15), label=rev(var.label) ); # Y-axis

abline(v=seq(-30,6,0.01), col="gray97", lty=1, lwd=1);

abline(h=c(2:4,6:8,10:12,14:16), lty=2, lwd=0.5, col="white");

#abline(v=seq(-5,5,0.5), col="white", lty=1, lwd=1);
abline(v=seq(-20,5,5), col="white", lty=1, lwd=2);


# 90% credible intervals
segments(par.mod2[7:10,2], c(15,11,7,3), par.mod2[7:10,3], c(15,11,7,3), lty=1, col="black", lwd=3);

# Means
points(x=par.mod2[7:10,1], y=c(15,11,7,3), pch=22, cex=1.2, col="black", bg="black");


